################################################################################
#################10.Table 10: Fixed effects with alt. income diff.##############
################################################################################

#####Load data from Data_Generation_CPS#####

pols <- readRDS('pols_bjps')

#####Load plm package#####
#install.packages('plm')
library(plm)
#install.packages('pcse')
library(pcse)
#install.packages('lmtest')
library(lmtest)
#install.packages('dplyr')
library(dplyr)

#####Alternative income differentiation Operationalizations#####

gcova <- plm(scale_econsdw ~ scale_dispgini + 
               scale_gcov_inc1 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_gcov_inc1 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_gcov_inc1: scale_socioeconomicsal +
               scale_gcov_inc1:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")
coeftest(gcova, vcov = function(x) vcovBK(x, cluster="group"))

gcovb <- plm(scale_econsdw ~ scale_dispgini + 
               scale_gcov_inc2 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_gcov_inc2 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_gcov_inc2: scale_socioeconomicsal +
               scale_gcov_inc2:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")
coeftest(gcovb, vcov = function(x) vcovBK(x, cluster="group"))

gcovc <- plm(scale_econsdw ~ scale_dispgini + 
               scale_gcov_inc3 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_gcov_inc3 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_gcov_inc3: scale_socioeconomicsal +
               scale_gcov_inc3:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")
coeftest(gcovc, vcov = function(x) vcovBK(x, cluster="group"))

gcovd <- plm(scale_econsdw ~ scale_dispgini + 
               scale_gcov_inc4 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_gcov_inc4 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_gcov_inc4: scale_socioeconomicsal +
               scale_gcov_inc4:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")
coeftest(gcovd, vcov = function(x) vcovBK(x, cluster="group"))

gcove <- plm(scale_econsdw ~ scale_dispgini + 
               scale_gcov_inc5 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_gcov_inc5 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_gcov_inc5: scale_socioeconomicsal +
               scale_gcov_inc5:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")
coeftest(gcove, vcov = function(x) vcovBK(x, cluster="group"))

#####GCOV Point Estimates#####

fu1 <- as.data.frame(summary(gcova, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) 

fu2 <- as.data.frame(summary(gcovb, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(gcovc, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(gcovd, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(gcove, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Income Dif.", "Gini*Salience", "Income Dif.*Salience", "Gini*Income Dif.*Saience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini","Income Dif.*Salience", "Gini*Income Dif.", "Gini*Salience","Gini*Income Dif.*Saience", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))

fit1e <- c(summary(gcova)$r.squared, summary(gcova)$fstatistic$statistic, summary(gcova)$fstatistic$p.value)
fit2e <- c(summary(gcovb)$r.squared, summary(gcovb)$fstatistic$statistic, summary(gcovb)$fstatistic$p.value)
fit3e <- c(summary(gcovc)$r.squared, summary(gcovc)$fstatistic$statistic, summary(gcovc)$fstatistic$p.value)
fit4e <- c(summary(gcovd)$r.squared, summary(gcovd)$fstatistic$statistic, summary(gcovd)$fstatistic$p.value)
fit5e <- c(summary(gcove)$r.squared, summary(gcove)$fstatistic$statistic, summary(gcove)$fstatistic$p.value)


fit.statse <- rbind(fit1e, fit2e, fit3e, fit4e, fit5e)

print(mean(fit.statse[,1])) #R-Squared
print(mean(fit.statse[,2])) #Adjusted R-Squred
print(mean(fit.statse[,3])) #F-Statistic
print(mean(fit.statse[,4])) #p-value of F-Statistic

#####Fixed Effect Theil#####

theila <- plm(scale_econsdw ~ scale_dispgini + 
                scale_gtheil_inc1 + 
                scale_socioeconomicsal +
                scale_dispgini : scale_gtheil_inc1 + 
                scale_dispgini*scale_socioeconomicsal +
                scale_gtheil_inc1: scale_socioeconomicsal +
                scale_gtheil_inc1:scale_socioeconomicsal : scale_dispgini +
                scale_galtansdw +
                scale_enep +
                scale_unemp,
              data = pols, 
              na.action = na.omit, 
              model = "within")

coeftest(theila, vcov = function(x) vcovBK(x, cluster="group"))

theilb <- plm(scale_econsdw ~ scale_dispgini + 
                scale_gtheil_inc2 + 
                scale_socioeconomicsal +
                scale_dispgini : scale_gtheil_inc2 + 
                scale_dispgini*scale_socioeconomicsal +
                scale_gtheil_inc2: scale_socioeconomicsal +
                scale_gtheil_inc2:scale_socioeconomicsal : scale_dispgini +
                scale_galtansdw +
                scale_enep +
                scale_unemp,
              data = pols, 
              na.action = na.omit, 
              model = "within")

coeftest(theilb, vcov = function(x) vcovBK(x, cluster="group"))

theilc <- plm(scale_econsdw ~ scale_dispgini + 
                scale_gtheil_inc3 + 
                scale_socioeconomicsal +
                scale_dispgini : scale_gtheil_inc3 + 
                scale_dispgini*scale_socioeconomicsal +
                scale_gtheil_inc3: scale_socioeconomicsal +
                scale_gtheil_inc3:scale_socioeconomicsal : scale_dispgini +
                scale_galtansdw +
                scale_enep +
                scale_unemp,
              data = pols, 
              na.action = na.omit, 
              model = "within")

coeftest(theilc, vcov = function(x) vcovBK(x, cluster="group"))

theild <- plm(scale_econsdw ~ scale_dispgini + 
                scale_gtheil_inc4 + 
                scale_socioeconomicsal +
                scale_dispgini : scale_gtheil_inc4 + 
                scale_dispgini*scale_socioeconomicsal +
                scale_gtheil_inc4: scale_socioeconomicsal +
                scale_gtheil_inc4:scale_socioeconomicsal : scale_dispgini +
                scale_galtansdw +
                scale_enep +
                scale_unemp,
              data = pols, 
              na.action = na.omit, 
              model = "within")

coeftest(theild, vcov = function(x) vcovBK(x, cluster="group"))

theile <- plm(scale_econsdw ~ scale_dispgini + 
                scale_gtheil_inc5 + 
                scale_socioeconomicsal +
                scale_dispgini : scale_gtheil_inc5 + 
                scale_dispgini*scale_socioeconomicsal +
                scale_gtheil_inc5: scale_socioeconomicsal +
                scale_gtheil_inc5:scale_socioeconomicsal : scale_dispgini +
                scale_galtansdw +
                scale_enep +
                scale_unemp,
              data = pols, 
              na.action = na.omit, 
              model = "within")

coeftest(theile, vcov = function(x) vcovBK(x, cluster="group"))

#####GTHEIL Point Estimates#####

fu1 <- as.data.frame(summary(theila, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) 

fu2 <- as.data.frame(summary(theilb, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(theilc, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(theild, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(theile, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Income Dif.", "Gini*Salience", "Income Dif.*Salience", "Gini*Income Dif.*Saience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini","Income Dif.*Salience", "Gini*Income Dif.", "Gini*Salience","Gini*Income Dif.*Saience", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))


rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))

#Model Fit Statistics#
fit1d <- c(summary(theila)$r.squared, summary(theila)$fstatistic$statistic, summary(theila)$fstatistic$p.value)
fit2d <- c(summary(theilb)$r.squared, summary(theilb)$fstatistic$statistic, summary(theilb)$fstatistic$p.value)
fit3d <- c(summary(theilc)$r.squared, summary(theilc)$fstatistic$statistic, summary(theilc)$fstatistic$p.value)
fit4d <- c(summary(theild)$r.squared, summary(theild)$fstatistic$statistic, summary(theild)$fstatistic$p.value)
fit5d <- c(summary(theile)$r.squared, summary(theile)$fstatistic$statistic, summary(theile)$fstatistic$p.value)


fit.statsd <- rbind(fit1d, fit2d, fit3d, fit4d, fit5d)

print(mean(fit.statsd[,1])) #R-Squared
print(mean(fit.statsd[,2])) #Adjusted R-Squred
print(mean(fit.statsd[,3])) #F-Statistic
print(mean(fit.statsd[,4])) #p-value of F-Statistic

#####Fixed Effect: SDMW#####

sdmwa <- plm(scale_econsdw ~ scale_dispgini + 
               scale_sdmw_inc1 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_sdmw_inc1 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_sdmw_inc1: scale_socioeconomicsal +
               scale_sdmw_inc1:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(sdmwa, vcov = function(x) vcovBK(x, cluster="group"))

sdmwb <- plm(scale_econsdw ~ scale_dispgini + 
               scale_sdmw_inc2 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_sdmw_inc2 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_sdmw_inc2: scale_socioeconomicsal +
               scale_sdmw_inc2:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(sdmwb, vcov = function(x) vcovBK(x, cluster="group"))

sdmwc <- plm(scale_econsdw ~ scale_dispgini + 
               scale_sdmw_inc3 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_sdmw_inc3 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_sdmw_inc3: scale_socioeconomicsal +
               scale_sdmw_inc3:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(sdmwc, vcov = function(x) vcovBK(x, cluster="group"))

sdmwd <- plm(scale_econsdw ~ scale_dispgini + 
               scale_sdmw_inc4 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_sdmw_inc4 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_sdmw_inc4: scale_socioeconomicsal +
               scale_sdmw_inc4:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(sdmwd, vcov = function(x) vcovBK(x, cluster="group"))

sdmwe <- plm(scale_econsdw ~ scale_dispgini + 
               scale_sdmw_inc5 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_sdmw_inc5 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_sdmw_inc5: scale_socioeconomicsal +
               scale_sdmw_inc5:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(sdmwe, vcov = function(x) vcovBK(x, cluster="group"))

#####SDMW Point Estimates#####

fu1 <- as.data.frame(summary(sdmwa, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) 

fu2 <- as.data.frame(summary(sdmwb, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(sdmwc, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(sdmwd, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(sdmwe, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Income Dif.", "Gini*Salience", "Income Dif.*Salience", "Gini*Income Dif.*Saience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini","Income Dif.*Salience", "Gini*Income Dif.", "Gini*Salience","Gini*Income Dif.*Saience", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))


rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))

#Model Fit Statistics#

fit1c <- c(summary(sdmwa)$r.squared, summary(sdmwa)$fstatistic$statistic, summary(sdmwa)$fstatistic$p.value)
fit2c <- c(summary(sdmwb)$r.squared, summary(sdmwb)$fstatistic$statistic, summary(sdmwb)$fstatistic$p.value)
fit3c <- c(summary(sdmwc)$r.squared, summary(sdmwc)$fstatistic$statistic, summary(sdmwc)$fstatistic$p.value)
fit4c <- c(summary(sdmwd)$r.squared, summary(sdmwd)$fstatistic$statistic, summary(sdmwd)$fstatistic$p.value)
fit5c <- c(summary(sdmwe)$r.squared, summary(sdmwe)$fstatistic$statistic, summary(sdmwe)$fstatistic$p.value)


fit.statsc <- rbind(fit1c, fit2c, fit3c, fit4c, fit5c)

print(mean(fit.statsc[,1])) #R-Squared
print(mean(fit.statsc[,2])) #Adjusted R-Squred
print(mean(fit.statsc[,3])) #F-Statistic
print(mean(fit.statsc[,4])) #p-value of F-Statistic
